Dynamics of an Ensemble of Noisy Bistable Elements with Global Time-Delayed 

Coupling 
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The dynamics of an ensemble of bistable elements with global time-delayed coupling under the 
influence of noise is studied analytically and numerically. Depending on the noise level the system 
undergoes ordering transitions and demonstrates multi-stability. That is, for a strong enough pos- 
itive feedback it exhibits a non-zero stationary mean field and a variety of stable oscillatory mean 
field states are accessible for positive and negative feedback. The regularity of the oscillatory states 
is maximal for a certain noise level, i.e., the system demonstrates coherence resonance. While away 
from the transition points the system dynamics is well described by a Gaussian approximation, near 
the bifurcation points a description in terms of a dichotomous theory is more adequate. 
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Noise induced hopping events in bi- or multi-stable 
systems form the basis of many interesting phenomena 
observed in physics, biology, chemistry, as well as social 
science. The study of such rate processes has thus been 
a subject of great interest and various techniques such 
as Langevin, Fokker-Planck and master equations have 
been used to describe the dynamics of stochastic systems 

In this article we generalize a well-studied stochastic 
model consisting of an ensemble of interconnected noise 
driven bistable oscillators by introducing time-delayed 
couplings. Such an extension is important since it has 
been realized that time-delays are ubiquitous in nature 
and often change fundamentally the dynamics of the sys- 
tem null. 

The dynamics of the network is studied numerically 
by using Langevin equations and analytically by two 
complementary mean field descriptions which are de- 
rived from the corresponding Fokker-Planck and master- 
equations, respectively. 

We assume that the bistable elements are highly in- 
terconnected, so that the connectivity can be approxi- 
mated by a global all-to-all coupling. In such a system a 
bistable element may for instance model the basic prop- 
erties of a neuron that can either be in a firing or a non- 
firing state, a person that can opt between two choices 
by means of individual voting, or a gene that is either 
expressed or non-expressed. The globally coupled sys- 
tem may then represent a highly interconnected neural 
network , a social group in which the individual voting 
behavior is influenced by opinion polls 0, or a genetic 
regulatory network , res pec tively. Other examples are 
catalytic surface reactions |lO| and allosteric enzymic re- 
actions im . 

The properties of globally coupled elements have been 
a subject of many studies 0, y, 0|- In particu- 
lar, Desai and Zwanzig 0] studied the synchronization of 
thermally activated bistable elements with instantaneous 
coupling and found an exact mean field solution in the 



thermodynamic limit N ^ oo, where N is the number of 
elements in the network. This system exhibits a second 
order phase transition to an ordered state with non-zero 
stationary mean field. The effect of a time-delayed cou- 
pling has been studied by Yeung and Strogatz 13] for a 
globally coupled network of periodic phase oscillators and 
Tsimring and Pikovsky 0| investigated the dynamics of 
a single bistable element driven by noise and time-delayed 
feedback. 

Here, as in it is assumed that the time delays 

between the bistable elements are uniform. Such an ap- 
proximation is justified in certain neural networks, whose 
time-delays are remarkably constant, as suggested by re- 
cent findings 'iS"]. Similarly, for certain regulatory ge- 
netic networks the response time lag is determined by a 
single time constant 16]. 

Our prototype system for the study of rate processes 
in extended systems consists of N Langevin equations, 
each describing the overdamped noise driven motion of a 
particle in a bistable potential V = — -|- a;^/4, whose 
symmetry is distorted by a global coupling to the time- 



delayed mean field X{t — T)=N ^ Y^^=i ~ ''')' 
x^{t) = x^it) - x,{tf + eX{t - t) 



(1) 



where r is the time delay, e is the coupling strength of 
the feedback and D denotes the variance of the Gaussian 
fiuctuations (,{t), which are mutually independent and 
uncorrelated {£,t{t)£,j{t')) = 5^j6{t ~ t'). 

We first study system ^ numerically. The simulations 
are carried out using a fixed-step fifth-order Runge-Kutta 
method with linear interpolation for the evaluations at 
intermediate steps required for the delayed variables. If 
not otherwise stated, the time-step and the number of 
elements is At = 0.01 — 0.05 and N = 2500, respectively. 

For £ = 0, the elements are decoupled from each other. 
They jump from one potential well to the other randomly 
and independently of each other. Therefore, in this case 
the mean field X = 0. For small ]e], the mean field 
remains zero. At a certain e — Est > ^ which depends 
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on the noise intensity D, but is independent of the time 
delay, the system undergoes a second order (continuous) 
phase transition and adopts a non-zero stationary mean 
field. 

For a negative feedback, a transition to a periodically 
oscillating mean field solution is observed at a certain e = 
Eosc- < 0. Here and for the rest of this paper a (— /+) 
index means that the corresponding value is associated 
with a negative/positive feedback. 

Above a certain noise level Dh the transition at Sosc- is 
second order as well. However, for D < Dh the system 
exhibits a first order (discontinuous) transition associ- 
ated with hysteretic behavior. The critical noise strength 
Dh depends on the time delay and is Dh = 0.07 for 
T = 100. 

For large time delays t ^ tk [tk is the inverse 
Kramers escape rate from one well into the other), de- 
pending on the initial state the system adopts one of the 
many accessible oscillatory states featuring different peri- 
ods. Even for a positive feedback, besides the stationary 
solution several oscillatory states with periods T < t are 
observed for e > £osc+ ^ ^st- If the feedback is negative, 
the system only has oscillatory non-trivial solutions. The 
observed periods are T < 2r for £ < £osc-- 

The simulations show that for a negative feedback all 
oscillating states have a vanishing time-averaged mean 
field {X)t- However, for a positive feedback besides sym- 
metric periodic solutions with {X)t — 0, states with sig- 
nificant non-zero temporal mean are possible. 

In order to theoretically study the dynamical proper- 
ties of a globally coupled set of noisy bistable elements 
(with no time delay), Desai and Zwanzig derived a 
hierarchy of equations for the cumulant moments of the 
distribution function from the multi-dimensional Fokker- 
Planck equation for the joint probability distribution for 
all elements. For large noise intensities, when the statis- 
tics of individual elements are approximately Gaussian, 
this hierarchy can be truncated. Applying this approach 
to our system yields the following set of equations for the 
mean field X and the variance M = N^^ ^ 

X{t) = X{t)-X^{t)-3X{t)M{t)+eX{t-T), 

1 (2) 
-M{t) = M{t)-3X^{t)M{t)-3M^{t)+D. 

To compare the predictions of the Gaussian approxi- 
mation (O with the original Langevin model we cal- 
culate the power of the main peak Ppeak in the power 
spectrum of X{t). It is proportional to the amplitude of 
the mean field oscillations and can thus be used to ana- 
lyze the Hopf bifurcation, which describes the transition 
to the oscillating mean-field states. The pitchfork bifur- 
cation to the non-zero stationary mean field for a large 
enough positive coupling is characterized by the depen- 
dence of the temporal mean of the mean field {X)t on 
the system parameters. 
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FIG. 1: The peak power Ppeak and {X)t as a function of the 
coupling strength for the Langevin model (crosses), the Gaus- 
sian approximation (dashed line), where the double line indi- 
cates hysteretic behavior, and the dichotomous theory (solid 
line). The constant noise strength is D=0.1 and the time de- 
lay is T = 100. For For X = and D = 0.1 the Kramers time 
is Tk ~ 54.1. 
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FIG. 2: Phase diagram for r = 100 of the Langevin model 
(crosses), the Gaussian approximation (dashed lines) and the 
dichotomous theory (solid lines and dotted lines). The solid 
line and the dotted line respectively depict the primary solu- 
tion and the higher order solutions of Eq. Q and Q . Phases 
separated by double lines indicate hysteretic behavior. The 
inset shows how the Hopf bifurcation line Eosc- (P) varies with 
the time delay r. For X = and D < 0.3 the Kramers time 
is Tk > 10. 

For T = 100 the peak power Ppeak and the temporal 
mean {X)t resulting from the Gaussian approximation 
and the Langevin model are shown in Fig. ^as a function 
of the coupling strength e. The phase diagrams of these 
models are shown in Fig. |21in the {D, £)-parameter plane. 

Fig. n shows that away from the bifurcation points 
the Gaussian approximation describes the Langevin dy- 
namics correctly. However, near the transition points 
the Langevin dynamics is strongly non-Gaussian even for 
large noise temperatures. For instance, the Gaussian ap- 
proximation predicts that both bifurcations are subcrit- 
ical for the entire noise range D = 0.03 — 1.0 considered 
in this study (see Fig. [2Jl, while in the original Langevin 
model the bifurcations are subcritical only for D < 0.07. 

Including higher-order cumulant equations only leads 
to a very slow convergence towards the true solution of 
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the Langevin model. Thus, in order to describe the be- 
havior of the system near the bifurcation points, we ap- 
ply a complementary dichotomous approximation, which 
is valid in the limit of small noise, where the character- 
istic Kramers transition time Tif ^ 1. In the dichoto- 
mous approximation intra-well fluctuations of Xi are ne- 
glected. Thus each bistable element can be replaced by a 
discrete two-state system which can only take the values 
xi^2 — ±1- Then the collective dynamics of the entire 
system is described by the master equations for the oc- 
cupation probabilities of these states ni^2- This approach 
has been successfully used in studies of stochastic reso- 
nance and coherence resonance [e.g. 0, ll3.JrR H^ . For 
instance, using this approach Jung et al. j3| found sta- 
tionary mean field solutions in a globally coupled, time 
delay free network of bistable elements. 

The dynamics of a single element is determined by the 
hopping rates pu and P2i, i.e., by the probabihties to 
hop over the potential barrier from xi to X2 and from 
X2 to xi, respectively. In a globally coupled system, ni,2 
and pi2,2i are identical for all elements. Then the master 
equations for the occupation probabilities read 



Til 



(3) 



In the dichotomous approximation the mean field X = 
Xini+X2n2 ~ n2 — ni, and making use of the probability 
conservation rai + n2 = 1 we obtain the equation for the 
mean field 



Xit) =Pl2 -P21 - {P21 +Pl2)X{t). 



(4) 



The hopping probabilities pi2,2i are given by the 
Kramers transition rate for the instantaneous poten- 
tial well, which in the limit of small noise D and coupling 
strength e reads [cf. Q| , 



Pi2,2i = -z exp 



27r 



41) 



(5) 

A linear stability analysis of Eq. Q near the trivial 
state X = yields the transcendent equation for the 
complex eigenvalue A, 



(6) 



For a positive coupling this equation always has a real 
solution. At a certain critical couphng est = 4D/ (4— 3£>) 
the eigenvalue becomes positive indicating the pitchfork 
bifurcation observed in the Langevin system . Besides 
this real solution, Eq. © possesses an infinite number 
of complex solutions corresponding to oscillating mean 
fields. However, only a finite number of them corresponds 
to unstable modes at finite r and e as we will see below. 
The critical values e of the corresponding instabilities are 



found by substituting A = /i -I- iw into Eq. ® , separating 
real and imaginary part and setting /i = 0: 



V2 



UJT — exp(— l/4_D)Ttanwr 

TT 

Est 



cos OJT 



(7) 
(8) 



This set of equations has a multiplicity of solutions, indi- 
cating that multi-stability occurs in the globally coupled 
system beyond a certain coupling strength. For finite 
time delays and positive coupling, besides the station- 
ary solution, several oscillatory states with periods 
close to but slightly larger than r/fc are observed for 
£ > ^osc-t- — li2,...}, where the transition points 
are ordered as follows, < Est < £osc+ < ^osc+ • • ■ If the 
feedback is negative, the system has oscillatory solutions 
with periods T; close to but slightly larger than 2t/ (2/+1) 
for £ < e\^^_ = 0, 1, . . .}, where > e^^c- > eL- ■ • ■ 

Let us now discuss the bifurcation properties in the 
limit of large and small time delays as well as vanishing 
noise and compare them with those of a single oscillator 
system. The critical coupling e^t of the pitchfork bifurca- 
tion is time delay independent and goes to zero for van- 
ishing noise. However, the critical coupling of the Hopf 
bifurcation depends on the time delay (see Fig|3 inset). 
As the time delay increases, the maximum of the primary 
Hopf bifurcation line e\^^_ approaches the origin in the 
(e, D) plane meaning that oscillations may occur at an 
arbitrary small feedback strength for the properly tuned 
noise level. This should be contrasted to the dynamics of 
a single noise-free oscillator with time-delayed feedback 
that only exhibits oscillations at strong negative feedback 
(e < —1). For very small time delays t ^ 0, the critical 
coupling strength e^^-^ =Foo. 

In order to compare the predictions of the dichotomous 
model with the Langevin dynamics, the peak power Ppoak 
and the temporal mean {X)f, resulting from the dichoto- 
mous theory, are also plotted in Fig. The phase dia- 
gram for the dichotomous theory is shown in Fig. |21 

Fig. n and 121 show that the dichotomous theory agrees 
with the Langevin dynamics quite well for small noise 
in the range D « 0.07 — 0.3 in the neighborhood of the 
bifurcation points. The theory also correctly describes 
the bifurcation type. Indeed, the dichotomous theory 
predicts accurately the noise strength Dh {— 0.07 for 
T = 100) at which the Hopf bifurcation changes from 
supercritical to subcritical. However, for very small D 
the Kramers time becomes very large, and the accuracy 
of numerics becomes insufficient for a comparison with 
the theory. 

Let us point out that the system studied in this paper 
exhibits the phenomena of coherence resonance [e.g. |20| 
and array-enhanced resonance. Since both Kramers ran- 
dom switching rate p (see Eq. and the frequency of the 
oscillatory states / = w/(27r) (see Eq. [7|) depend on the 
noise strength, i.e., p — p{D) and / = f{D), the noise 
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FIG. 3: The normalized coherence of the oscillatory states 
/3 for systems of different size A'^ with r — 100. (a) (3 as 
a function of the noise strength D for e = 0.2. (b) /3 as a 
function of the coupling strength e for D = 0.1. For £ < 
and e > the spectral peak frequency is fp ~ 0.5 1/r and 
fp ~ 1.0 1/r, respectively. The dotted vertical line depicts 



can be tuned so that the random hopping between the po- 
tential wells of the bistable oscillators synchronizes with 
the periodic modulation of the mean field. This statisti- 
cal synchronization takes place when f — p/2 p^ . where 
the regularity of the oscillatory state becomes maximal. 
Here, this regularity is quantified by /3 = Hfp/Af, where 
H is the hight of the spectral peak at fp and A/ is the 
half width of the peak. The coherence measure /? as a 
function of the noise strength is shown in Fig. El^a). We 
observe that the regularity of the oscillatory states in- 
creases with increasing TV, a property which was reported 
for other systems and is sometimes referred to as array- 
enhanced resonance [2^ . Interestingly, the enhancement 
of the temporal regularity with increasing system size 
is only observed for macroscopic mean field oscillations, 
while the inverse holds for "subcritical coherence" . That 
is, the coherence observed in the power spectra of subcrit- 
ical mean field fluctuations (i.e., for |e| < |eosc±|) decays 
inversely proportional to the number of elements in the 
network, and becomes negligible for N > 10. This is 
shown in Fig. Olb). Qualitatively, the same dependency 
on the system size is found if the delayed average does 
not include the delayed element itself, i.e., the element 



N-l 



Xi is coupled to Xi{t — t) = ^ 

In summary, we have shown that a network of noisy 
bistable elements with global time-delayed coupling pos- 
sesses a multiplicity of stable oscillatory states for both 
positive and negative feedback in addition to a non-zero 
stationary mean field for a strong enough positive feed- 
back which also occurs in a non-delayed system. These 
novel oscillatory states have a maximum regularity for 
a certain noise strength. The bifurcations of the trivial 
equilibrium are well described by the dichotomous the- 
ory in the limit of small noise and coupling strength. Far 
away from the bifurcation points the mean-field prop- 
erties of the system are well described by the Gaussian 
approximation. However, the quantitative theory for the 
large noise strength near the bifurcation points is still 
lacking. In this paper the effect of uniform time de- 



lays on the dynamics of a globally coupled network of 
bistable elements has been studied. However, many real 
networks have sparse coupling and non-uniform time de- 
lays. These properties should thus be included in a more 
general description. Preliminary results of our simula- 
tions with nonuniform time-delays suggest, that the bi- 
furcation properties do qualitatively not change for a 
wide range of Gaussian distributed time delays, which 
substantiates the generic nature of the here considered 
model. This issue will be discussed in detail in an up- 
coming article. 
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